This is version 0.17.04.05.

Background

Fill in some background info.

Requirements

All analyses require the R software (v3.3) for data simulation, processing, and summarizing model results, and the Stan software for Hamiltonian Monte Carlo (HMC) simulation. Please note that some of the R code below may not work with older versions of JAGS due to some changes in the ways that arrays are handled.

We also need a few packages that are not included with the base installation of R, so we begin by installing them (if necessary) and then loading them.

if(!require("tvvarss")) {
  if(!require("devtools")) {
    install.packages("devtools")
    library("devtools")
  }
  devtools::install_github("eric-ward/tvvarss")
  library("tvvarss")
}

Step 1: Simulate data

The tvvarss package has a function simTVVAR to simulate the process component of a TVVARSS model (i.e., it does not add observation error). We use that to simulate multiple realizations of a given process and then see how well the true parameters can be estimated.

Food web topologies

The food web topology in a TVVARSS model is expressed via the matrix \(\mathbf{B}_t\). Specifically, interactions are expressed as the effect of column on row; the diagonals indicate the strength of density-dependence. All elements of the matrix corresponding to no interaction are set to 0.

tvvarss is designed to work with symbolic representations within \(\mathbf{B}_t\), based on the following codes:

  • "dd": density-dependence (this is implied in TVVARSS models)
  • "td": top-down
  • "bu": bottom-up
  • "cf": competitive/facilitative

We show four different examples of simulated food web topologies.

Ex 1: Linear food chain

Here there are 4 tropic levels stacked in a linear food chain from primary producers PP at the bottom to tertiay consumers TC at the top.

lvls <- c("TC","SC","PC","PP")
cat(paste0(paste0(lvls[1:3],"\n|\n",collapse = ""),lvls[4]))
## TC
## |
## SC
## |
## PC
## |
## PP

Here is how to express that topology in a matrix form that tvvarss will understand.

nn <- length(lvls)
## initial conditions for B_t
B0_lfc <- matrix(list(0),nn,nn)
dimnames(B0_lfc) <- list(lvls,lvls)
## diagonal elements = density-dependence
diag(B0_lfc) <- rep("dd",4)
for(i in 1:(nn-1)) {
  B0_lfc[i,i+1] <- "bu"
  B0_lfc[i+1,i] <- "td"
}
## inspect B0
B0_lfc
##    TC   SC   PC   PP  
## TC "dd" "bu" 0    0   
## SC "td" "dd" "bu" 0   
## PC 0    "td" "dd" "bu"
## PP 0    0    "td" "dd"

We can now simulate a TVVAR process based on this expression of the food web topology. In addition to the matrix specifying the topology, the simulator needs to know the length of time series and some information about the variances of the process errors for \(\mathbf{B}_t\) and the states \(\mathbf{x}_t\).

Here is a 30-unit long TVVAR model.

TT <- 35
## simulate
set.seed(666)
lfc <- simTVVAR(B0 = B0_lfc, TT = TT,
                var_QX = seq(4)/40, cov_QX = 0,
                var_QB = 0.05, cov_QB = 0)

And we can plot the states over time.

par(mai=c(0.9,0.9,0,0), omi=c(0.1,0.1,0.1,1.5))
clr <- c("purple","darkred","blue","darkgreen")
matplot(t(lfc$states), type="l", lty="solid", lwd=2, xpd=NA,
        col=clr)
legend("right", legend=lvls, lty="solid", lwd=2,
       col=clr, inset=-0.2, xpd=NA)

Let’s fit a bunch and inspect them.

ee <- 10
ex1 <- vector("list",ee)
for(i in 1:ee) {
  ex1[[i]] <- simTVVAR(B0 = B0_lfc, TT = TT,
                       var_QX = seq(4)/40, cov_QX = 0,
                       var_QB = 0.05, cov_QB = 0)
  par(mai=c(0.6,0.9,0.3,0), omi=c(0.1,0.1,0.1,1.5))
  clr <- c("purple","darkred","blue","darkgreen")
  matplot(t(ex1[[i]]$states), type="l", lty="solid", lwd=2, xpd=NA,
          col=clr, main=paste0("Sim ",i), ylab="Log density")
  legend("right", legend=lvls, lty="solid", lwd=2,
         col=clr, inset=-0.2, xpd=NA)
}

OK, clearly some of them exhibit unrealistic, explosive population dynamics. Let’s develop a screen to toss those out.

## number of simulations
ee <- 10
## min log-density threshold
dens_min <- -6
## max log-density threshold
dens_max <- 10
ex1 <- vector("list",ee)
for(i in 1:ee) {
  tmp <- list(states=2*rep(dens_max,2))
  while(max(tmp$states) > dens_max | min(tmp$states) < dens_min) {
    tmp <- simTVVAR(B0 = B0_lfc, TT = TT,
                    var_QX = seq(4)/40, cov_QX = 0,
                    var_QB = 0.05, cov_QB = 0)
  }
  ex1[[i]] <- tmp
  par(mai=c(0.6,0.9,0.3,0), omi=c(0.1,0.1,0.1,1.5))
  clr <- c("purple","darkred","blue","darkgreen")
  matplot(t(ex1[[i]]$states), type="l", lty="solid", lwd=2, xpd=NA,
          col=clr, main=paste0("Sim ",i), ylab="Log density")
  legend("right", legend=lvls, lty="solid", lwd=2,
         col=clr, inset=-0.2, xpd=NA)
}